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Abstract: Contour maps are widely used to display estimates of spatial 
fields. Instead of showing the estimated held, a contour map only shows 
a fixed number of contour lines for different levels. However, despite the 
ubiquitous use of these maps, the uncertainty associated with them has 
been given a surprisingly small amount of attention. We derive measures 
of the statistical uncertainty, or quality, of contour maps, and use these 
to decide an appropriate number of contour lines, that relates to the un¬ 
certainty in the estimated spatial held. For practical use in geostatistics 
and medical imaging, computational methods are constructed, that can 
be applied to Gaussian Markov random helds, and in particular be used 
in combination with integrated nested Laplace approximations for latent 
Gaussian models. The methods are demonstrated on simulated data and 
an application to temperature estimation is presented. 

Key words: Excursions; Kriging; Latent Gaussian processes; Markov 
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1 Introduction 


Contour maps are often used in environmental statistics and applied spatial statis¬ 
tics to display estimates of continuous surfaces. One such example can be seen in 
Figure [l] where an estimate of the mean summer temperature in 1997 is shown for 
the US. 

One of the earliest documented uses of contour maps is from 1548, when water 
depth was displayed using contours in a map |Thrower, 2008 . Since then, drawing 


contour lines (sometimes also called isolines) has been a frequently used method 
for displaying 3D surfaces in print. In cartography, the number of contours, or 
the distance between them, is often chosen depending on the intended purpose of 
the map. The more contours that are drawn the greater detail can be extracted 
from the map, but the more information is also needed when drawing the map. 
Thus, the uncertainty in the information about the topography limits the number 
of contours that can be drawn. When contour maps are used to display estimated 
surfaces in modern spatial statistics, the number of contours used are meant to 
reflect the uncertainty in the estimate. Intuitively, one should be allowed to draw 
many contours if the uncertainty of the estimated surface is low and fewer contours 
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if the uncertainty is high. However, it has not been clear how the number of 
contours should be decided nor at what levels the contours should be drawn in 
order to accurately reflect the variability in the estimated surface. 

Even though contour maps are widely used, relatively little research has focused 
on quantifying their statistical properties. This problem is clearly related to the 
problem of finding uncertainty regions for contour curves, which was studied by 
Lindgren and Rychlik |1995 using theory of level crossings in Gaussian processes 


and later by Wameling and Saborowski |2003 using conditional simulations. Both 
these methods construct the uncertainty regions locally, without any control over 
the simultaneous coverage probabilities for the regions. This problem was recently 


addressed by French 2014 , using a simulation-based approach, and by Bolin and 


Lindgren 2015| using a unified approach for both the contour uncertainty prob¬ 


lem and the related problem of finding excursion sets of random fields. To our 
knowledge, the problem of deciding how many contour curves one should draw in 
a contour map has previously only been studied by Polfeldt |l999] , who proposed 
a method which we will come back to in Section 12.11 

The main contributions of this work are threefold. Firstly, we discuss different 
ways in which the statistical properties of contour maps can be understood. Based 
on these interpretations, we propose statistical quality measures that can be used 
to evaluate how appropriate a given contour map is for a random held. Secondly, 
we propose a strategy for deciding the number of contour levels to use in a contour 
map, and discuss how the actual contour levels should be chosen. Finally, we 
discuss the problem of discrete and continuous domain interpretations for contour 
maps and propose practical methods for interpolating discrete-domain calculations 
to continuous domains. This is an important topic, because the domain of interest 
for the analysis often is a continuous subset of M 2 , such as the unit square. The 
statistical computations have to be carried out for a finite number of locations, 
such as those on a regular lattice covering the domain. Thus, the contour map has 
to be constructed using some form of interpolation of discrete-domain calculations. 
The problem of analyzing a single contour curve for a surface obtained by linear 
interpolation of a kriging prediction on a lattice was briefly discussed by Wameling] 
120021 and the interpolation problem has to our knowledge not been studied at all 
for general contour maps before. 

Most of the theory that is presented later is applicable to any random held 
model. However, we need to restrict ourselves if we also want a computationally 
efficient method for applying the theory. To that end, we will focus on the class of 
latent Gaussian models, which is a popular model class that includes many of the 
standard models in spatial statistics Rue et ah, 2009]. A latent Gaussian model 


is specified hierarchically using a possibly non-Gaussian likelihood for the data 
y conditionally on a latent Gaussian held x(-). The model parameters, such as 
correlation ranges and measurement noise variances, typically have prior distribu¬ 
tions in a Bayesian setting, but whether or not priors are used does not affect the 
methods we will discuss here. The quantity of interest is the distribution of the 
latent held x(-) given the data, which typically is summarized using the posterior 
mean E[x(s)|y] as a point estimate (often called the kriging predictor in geostatis¬ 
tics) and the posterior variances V[x(s)|y] as a measure of uncertainty, or using a 
contour map of the posterior mean. 

The structure of the paper is as follows. In Section[2j we introduce some needed 
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Figure 1: Example of a contour map for the mean summer temperature in 1997 
in the US. The temperature changes due to topography in the western part is 
apparent. The challenge is to quantify our level of trust in the level of detail in 
the contour curves. 


notation and give a precise definition of contour maps. The section then contains 
three subsections that contain the main results of the paper: Section 2T] introduces 
the method by Polfeldt [l999| and generalizes it to a a new method for quantifying 
and displaying the uncertainty in a contour map. Section 2.2 introduces three 
different quality measures for contour maps, based on three different statistical 
interpretations of contour maps. Finally, Section T3 discusses how the contour 
levels can be chosen. Section[3]presents the practical aspects of how to compute the 
quality measures for discrete domains and presents a method for interpolating the 
discrete domain calculations to a continuous domain. Section |4] demonstrates the 
methods using simulated data and Section [5] shows an application to temperature 
estimation. Finally, Section [6] concludes with a discussion. All methods presented 


and Lindgren, 
Becker et ah 

2016j. The outline of the US was obtained with the maps package 
2016 and the Bayesian inference was performed using the INLA 

package | Lindgren and Rue 

2015, 

Lindgren et al., 

2011, 

Rue et al.[ 

2009, 

Martins 

et al., 

2013 . 


2 Contour maps for random fields 

In order to quantify uncertainty of contour maps for random fields, we first need 
to more clearly define what such a contour map means. Only then can questions 
about formal uncertainty quantification and appropriate numbers of contour levels 
be considered. One of the most important points is that in presenting results from 
spatial statistical inference, contour maps are based on spatial point estimates 
/(•), or functionals, of properties of the random field distribution, such as the 
posterior mean E[x(-)|y], and not on realizations of the held itself. The contour 
map is thus also itself a functional of the distribution of the held. There is no 
single choice of functional /(■) to base the contour map on, and the choice affects 
the interpretation of both the contour map itself, and any associated uncertainty 
measure for it. 
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Let 12 be the domain of interest for the analysis. The results will depend on 
joint properties of the stochastic process across the entire domain, so the choice of 
12 is important for the interpretation of the results. In order to avoid unnecessary 
theoretical complications, we assume that f2 is an open set with a well-defined area 
|Q < oo. If the domain of interest is not open, we take the interior and ignore the 
behavior on the boundary. 

First, we revisit some notation from Bolin and Lindgren |2015] for excursion 
sets and contour sets of fixed functions / (s), s 6 fi, This notation is also applicable 
to the true contours of realizations of the random held x(s), s G f2. 


Definition 2.1 (Excursion and contour sets for functions). Given a function f : 
12 i —y K, the positive excursion set for a level u is defined by Af(f) = {s 6 12; /( s) > 
u}. The negative excursion set is similarly defined by A~(f) = {s G 12; /( s) < u}. 
The contour curve A^ for a level u is given by A^(f) = (Af(f)° U A~(f)°) c , where 
B° denotes the interior of the set B and B c denotes the complement. 


It is important to note that we do not define the contour curve (which is a 
set in general) A c u as the set of locations {s : /( s) = u}, but rather as the set of 
all continuous and discontinuous level u crossings. This is different from how, for 
example French |2014 defines contour curves, and the reason for using this slightly 
more complicated definition is that we want the theory to be applicable also for 
discontinuous fields. Such fields occur frequently in environmental statistics, for 
example as soon as we have a discontinuous covariate for the mean value of the held. 
The contour set definition is applicable even to realizations of exotic random helds 
such as completely independent Gaussian noise processes, where the contour set for 
any finite level u becomes the entire domain, 12. Contours for such random helds 
are not very useful, so without any practical limitation we can restrict the analysis 
to functions and realizations of random helds that are piecewise continuous, i.e. the 
probability that the set of discontinuity points is a null set is 1. Also noteworthy, 
but producing no practical problems, is that in the presence of discontinuities, 
contour sets for different levels are not necessarily distinct. 

We can now give a precise definition for contour maps of a function /(•). 


Definition 2.2 (Contour maps). Let f : 12 i-> M be a function. A contour map 
with K contour levels u\ < U 2 < ... < uk is defined as the collection of contour 
curves A^f),... ,A^ K (f) and associated level sets Gk = {s : Uk < f(s) < Uk+i}, 
for 0 < k < K, where we define w 0 = — oo and uk+i = oo. We denote this contour 
map by Cf(ui,... ,uk )• 


An illustration of contour curves and the associated contour map for a specihc 
function /(•) is given in Figure [2j showing the level sets Gk- Further details on 
how the different definitions affect the interpretation of the contour sets and maps 


on continuous domains are given in Section 3.2 


We now turn to the problem of interpreting contour maps of functions /(■) 
as information about a random held x(-). Most commonly, a contour map for a 
random held x(-) is constructed by taking /( s) to be the posterior mean E[x(s)|y] 
of the held, but other contour maps can be obtained from the posterior median, or 
some other hxed quantity. The generated contour map does not by itself contain 
information about the variability in x(s)|?/. Therefore, we will now focus on ways 
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Figure 2: A function /(•) and its contour curves for levels Uk- 1 , Uk, Uk+i and Uk +2 
(left), and the associated contour map Cf(ui,... ,Uk ) (right), with the level sets 
Gk- 1 , Gk, and Gk+i labeled. In the right panel, the color for each region Gk is 
chosen as the color corresponding to the level u% = (uk +1 — Uk)/ 2. 


to measure how appropriate a given contour map is for a given random held. Even 
though we primarily have the situation of looking at a contour map for a posterior 
distribution in mind, we will simplify the notation by dropping the dependence on 
the data and consider the more general problem of quantifying how appropriate a 
contour map Cf(ui,..., uk) is for some general (typically non-stationary) random 
held x(-). 


2.1 Probabilistic contour map functions 


An intuitively natural interpretation of a contour map for a random held x(-) is that 
x(s) should stay between and Uk+i for all s G Gk, for every k — 0,..., K. With 
the interpretation, the joint probability P(w fc < x(s) < Uk+ i,s G Gk, 1 < k < K) 
should be large if the contour map is appropriate for the random held. It is there¬ 
fore tempting to take this joint probability as a quality measure. Unfortunately, 
this probability will be low, or even zero in most cases, since the marginal prob¬ 
abilities P (uk < x(s) < Uk+ 1 ) for any s close to the boundary of the set G k will 
be small. Hence, this interpretation is too strong in general and we need differ¬ 
ent interpretations of the contour map in order to construct meaningful quality 


measures. 


In order to seek a practical method for determining an appropriate contour 
level spacing, Polfeldt [1999 proposed a method based on analyzing the marginal 
probabilities 

p( s) = P {u k < x(s) < u k+ i), for k s.t. s G G k , (1) 

for a Gaussian random held with mean value function /i(s) and standard devia¬ 
tions cr(s). In our notation, Polfeldt’s analysis was based on the fact that if the 
distribution for x is Gaussian, these probabilities can be written as 


p( s) = $ 


(1 - r(s))(u fc+ i - u k ) 


as 




r(s)(ufc+i - u k ) 


as 


( 2 ) 


where $ denotes the standard Gaussian distribution function and 

Ms) 


r s = 




The ratio (uk+i — Uk)/cr( s) is generally dependent on the spatial location s. 
However, one can compute a common ratio q = (uk+i — Uk)/cr 0 for the entire 
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contour map if cr(s) can be approximated by a constant <7o and if the contour 
levels Uk are equally spaced (in the sense that Uk+i — Uk is independent of k). If 
cr(s) is replaced by cro also in (|2j) , one can plot p as a function of r for different 
values of q. Because r measures how close one is to a level contour, p( s) will 
always be close to 0.5 for r close to zero, but it will quickly increase towards 1 if 
q is chosen high enough. 

Polfeldt argued that a;(s) falls within the contour bands for most values of s if 
p(r) is close to one for most values of r, and one should therefore choose the number 
of contours so that q is large enough for this to occur. Because q is determined 
by the spacing between the contours and the average standard deviation of the 
kriging estimator, this is a reasonable conclusion. 

However, there are many approximations in Polfeldt’s argument: the strat¬ 
egy only works for Gaussian posterior distributions (although it is straightforward 
to extend it to other distributions), it assumes that the kriging errors are well- 
approximated by a constant, and it does not take the spatial dependency of the 
data into account. In practical applications, these approximations are often diffi¬ 
cult to justify, and we will therefore now focus on extending this method to relax 
(actually remove) the assumptions. The idea we will use is to extend the excursion 
functions introduced by Bolin and Lindgren 2015 to contour maps. 


In order to extend the method, we first need some notation. The following 
definition of the contour-avoiding set for a random held is taken from Bohn and 
Lindgren 12015]. 


Definition 2.3 (Contour-avoiding set). Let x( s), s E LI be a random field, let 1 — a 
be a predefined error probability, and let u be the value for a contour curve. The 
1 — a contour-avoiding set is defined as the union E ua (x) = M+ a (x) U M~ a (x). 
Here, 


( M u,a( x ), M u,a( x )) = argmax{| D UD + \ : P (D C A u (x), D + C Affx)) > 1-a}, 

(D+,D~) 


is the pair of joint contour u excursion sets, where the sets (D + ,D ) are open and 
\D\ denotes the area of the set D. 


Recall that D~ C A~{x) means that x(s) < u for s £ D ~, and that D + C 
Affx) means that a;(s) > u for s £ D + . Thus, the pair of contour-avoiding sets is 
the largest pair of sets where, with probability 1 — a, one has x(s) > u in one set 
and x(s) < u in the other. 

A credible region for the level u contour curve is defined as E° a (x) = (E ua (x)) c , 
which is the smallest set such that with probability at least 1 — a, all level u 
crossings of x are in the set. Thus, this credible region has a well-defined global 
interpretation, and one should note that the definition is different from some other 
definitions for contour uncertainty regions in the literature, e.g. that of Lindgren 
and Rychlik 11995 . Now, because we are interested in joint properties of contour 


maps, we extend this definition to multiple contours as follows. 


Definition 2.4. Letx(s), s £ be a random field, let u = [u\,... ,uk) where u\ < 
U 2 < ... < Uk are some predefined levels, and let 1 — a be a given probability. The 
collection of joint contour-avoiding sets for the levels u, = (M‘ .... M*J, 
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where the sets Dk are disjoint and open and Ak = A Uk+i (X) fl Af k (X) = {s : Uk < 
x(s) < |_i } . The joint u contour-avoiding set is given by the union of these sets, 


C u , a (x) = [J M^ c 


Recall that Dk C Gk(x) means that Uk < x(s) < Uk+i for s e Dk■ Thus, the set 
< 7 ^ 0 ( 0 :) is the largest set so that with probability at least 1 — a, the random held 
x satishes the natural interpretation of the contour map for the locations in the 
set. Hence, we would like to choose the number of contour lines so that C UiQ ,(a;) is 
large. This procedure is therefore a natural extension of the Polfeldt procedure. 

Now, we can go further and calculate a variant of the excursion functions 
introduced by Bolin and Linclgren |2015j. The contour map function, 


F u (s) = sup{1 - a; s e C U)Ce j, 


is obtained for each s by identifying the smallest a such that s is contained in the 
contour-avoiding set C UjQ . The function P u (s) takes values between zero and one 
and can be used to visualize the ct-indexed family of functions C u a . Note that each 
C u a can be retrieved as the 1 — a excursion set of the function F u (s). Visualizing 
this function is again a natural extension of Polfeldt’s procedure of visualizing p( s) 
as a function of r, because we want to choose the number of contours such that 
F u ( s) increases rapidly when we move away from the contour lines. 


2.2 Global quality measures 

We will call a function P(x,Cf ) a contour map quality measure if it takes values in 
[0,1] and if P « 1 indicates that Cf, in some sense, is appropriate as a description 
of the distribution of the random field x, whereas F « 0 indicates that Cf is 
inappropriate. Given the contour map function, we define a simple contour map 
quality measure, Pq, as the normalized integral of the contour map function, 

P 0 (x, C f ) — --yj- F u ( s) ds. (3) 

Loosely speaking, this measure tells us the percentage of the total area for which 
the statement of the contour map holds. 

By revisiting which sets are involved in the joint probability calculations, we 
will construct two additional quality measures, with geometrically interpretable 
properties. 

2.2.1 The Pi quality measure 

One possibility is to require that the unions Gk-iAGk with high probability should 
contain all level Uk crossings of the process x. We define the contour map measure 
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Pi as the probability for this occurring. To ensure that all values of x(-) are 
covered by crossing levels, and for ease of notation, we define Uj = — oo for j < 0 
and Uj = oo for j > K. Because U = A~ (f) and U fL k Gi = A+ (/), we can 
write the measure as 

Pi(x,C f ) = P ( p){x(s) <u k , se A~ k i (f)} n (x(s) > U k , se A+ +i (/)} 

\k= 0 

Thus, the interpretation is that if P\ is close to one, then for any k — 1,..., K the 
probability is small for x taking the value u k outside the two regions G k _\ and G k 
bordering the level u k contour curve. The following equivalent formulation of the 
measure is useful for practical calculations, 

( K 

Pi(x, Cf) = P I P|{wfc_i < z(s) < u k+2 ,s G G k } 

\k= 0 

The asymmetry in that the lower bounds are u k - 1 while the upper bounds are 
u k+2 is caused by the simple fact that each G k lies between u k and u k+l . More 
specifically, s G G k -1 U G k is equivalent to that x(s) < u k for s G G k ~2 and 
that x(s) > u k for s G G k+ \. This means that x(s) < u k + 2 for s G Gt and that 
x(s) > for s G G k , which we can write as u k _ 1 < x(s) < u k+2 for s G G k . 
From this formulation, we see if P\ is large we have, with high probability, 
< x(s) < u k+ 2 for s G G k , i.e. that the locations in G k might even have true 
values as high as the next level up or as low as the next level down. Thus, this 
measure is based on a weaker restriction on how the held rc(s) can vary compared 
to the natural deterministic interpretation that requires the function to lie strictly 
between two adjacent levels. For the choice /(s) = E[a:(s)], Figure ^middle) shows 
five realizations of the level u k contour curve for #(•), and they lie mostly inside 
G k - 1 U G k , as desired. 

2.2.2 The P 2 quality measure 

The P\ focuses in the contour curves, resulting in joint probabilities for overlapping 
sets UG k for k — 1,... K. Another possibility is to focus on the disjoint areas 
G k for k = 0,..., K. Because G k contains all locations associated with values 
between u k and u k+ 1 , it is natural to associate the area G k with the midpoint level 
u l = ( U k+Uk+ 1)/2. To ensure a sensible interpretation for G 0 and Gk, if K > 1 we 
set uq = 2u\ — u 2 and uk +1 = 2 uk — uk-i when defining Uq and u e K , and if K = 1 
we set u 0 = u 1 - sup sen /(s) + inf se n /(s) and u 2 = u x + sup sen /(s) - inf sef2 /(s). 
The values are also used in general when deciding the plotting color for each 
G k set, as in Figure [I] and Figure fright). The idea is then that each G k with 
high probability should contain all level u e k crossings, and we define the contour 
map measure P 2 by the probability for this occurring, 

P 2 (x, C f ) = P (p|{x(s) <«i,sG^(/)}n{i(s) >u e k ,s G A+ +1 (/)} 

\k= 0 

This can be interpreted as the simultaneous probability for the level crossings of 
(ul, ..., u e K ) all falling within their respective level sets (G\, ..., Gk), which is 
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Figure 3: An expectation function /(s) = E[x(s)], with contour curves for levels 
Uk, u k-n an d u t+ 1 mar ked (left), five realisations of the Uk contour of x(-) 
(middle), and five realizations of the u k contour of x(-) (right). In the middle 
panel, most Uk contours are in the set Gk U G k +\, with joint probability for all k 
quantihed by the P\ measure. In the right panel, most u k contours are in the set 
Gk, with joint probability for all k quantified by the P 2 measure. 


then a simultaneous credible region collection with credibility probability P 2 {x, Cf). 
This is close to how one intuitively interprets contour maps. For the choice 
/(s) = E[x(s)], Figure fright) shows five realizations of the level u k contour curve 
for x(-), and they lie mostly inside Gk, as desired. 

If P 2 is large, then for any k = 0,..., K, the probability is low for x taking 
the value u% outside the region Gk that is associated with that value. We can also 
reformulate this measure to an expression more suitable for practical computation, 

P 2 (x, Cf) = P ( flK,.! < x(s) < u e k+1 , s £ Gk} 

\k =0 

where we interpret Uj = — oo for j < 0 and Uj = oo for j > K. From this 
formulation, we see that, with high probability, u k _ x < x(s) < u k+l for s £ Gk- 
Because Uk-i < u k _ x and u k+1 < Uk+ 2 , w e have that P 2 puts a stronger restriction 
on how the field can vary compared with the Pi measure. 

All three quality measures Pq, Pi, and P 2 take values between 0 and 1 and can 
be used to quantify how appropriate a given contour map is for a given stochastic 
field. The Po measure is linked to the spatially interpretable contour map function 
P u (•), whereas the P 2 measure gives the joint credibility probability for crossings 
of levels in between the ones displayed in a standard contour map. The measure 
Pi is similar in spirit to P 2l but is more permissive. 

2.3 Choosing the contour levels 

The quality measures presented above can be used to choose an appropriate num¬ 
ber of contours. This is done by first deciding on a credibility level, such as 90%, 
and then choosing the largest number of levels, K, so that the quality measure for 
the contour map is above the chosen credibility level. 

Given K one must also choose the contour levels Ui ,..., Uk- When doing this 
one can either permit arbitrary levels, or mandate a common level spacing, i.e. 
that the difference Uk+i — Uk is the same for all k. Given a restriction of using 
exactly K levels within the range of f(-), one choice, here referred to as Standard, 
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is to place the contour levels Uk with an even spacing between the minimal and 
maximal value of the function /. 

2014), 


The Standard choice is the default in contour in MATLAB MATLAB 


whereas the contour function in R |R Core Team 2013 by default uses a more 
unpredictable heuristic in order to ensure aesthetically pleasing level combinations. 
More specifically, for a given K, the method finds a sequence of about K+l equally 
spaced values which cover the range of /. The values are chosen so that they are 
1, 2 or 5 times a power of 10. We refer to this choice as Pretty, and examples of 
both Pretty maps and Standard maps are shown for simulated data in Section [4j 

For Standard contour maps, there is a direct relationship between K and the 
level spacing. Because the contour levels can cover more than the range of /, this 
is not the case for Pretty contour maps. Therefore, the level spacing is of more 
direct interest than K when comparing different contour maps. See Figure [7] for 
examples of pretty contour maps with the same K but with different level spacing. 

If the restriction of a common level spacing were to be removed, an alterna¬ 
tive would be to choose the non-equally spaced levels that maximize the quality 
measures subject to the restriction that the levels cover the range of /. Compu¬ 
tationally, a cheap approximate solution to this can be obtained using the bounds 
for the P measures presented in Section |3.1[ However, due to the lack of a clear 


interpretation of the resulting contour maps, we do not pursue this further. 


3 Practical considerations and continuous domains 


We have now presented several quality measures for contour maps, including an 


extension to the Polfeldt 1999 procedure for choosing the number of contour lines. 


However, these tools would not be of much use unless we also have a method for 
computing them. Thus, in this section we go through the practical details of how 
to calculate the different quality measures. 

A common scenario for when contour maps are used is when the domain H is 
a continuous subset of M 2 , such as the unit-square. In order to do any calculations 
in this case, one has to discretize the problem. We therefore start with the simpler 
problem of when H is discrete in Section 3.1| After this, we briefly describe how 
one can approximate the continuous-domain case using corrections to the discrete- 


domain computations in Section 3.2 


3.1 Calculating the measures for discrete domains 

Assume that the process is defined on a discrete domain with n locations {sx,..., s n }, 
such as a regular lattice, so that x = {x ±,... ,x n }. Furthermore, let 7r(x) be the 
joint probability density function for the process, and let tt (xi) denote the marginal 
density function for location i. 

Simple upper bounds for the measures Pi and P 2 can be obtained using the 
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marginal probabilities for the X{ variables. Let 


Pi = 



r u j+ 2 


= P{Uj -1 

<C <C Zlj-^2') 

/ 7T( 

[ Xi)dx, 




Juj- 1 





r u i+ 1 


p\ -- 

= PK -1 

<Xi< u e j+1 ) = 

/ 7T( 

[xi)dx, 




Ju j - 1 



where j is the index such that .s t £ Gj. The joint measures now satisfy Pi < 
min(p-) and P 2 < min(p^). If we are interested in a contour map where P > 1 — a, 
we can use these bounds to quickly reject all contour maps where min(p^) < 1 — a 
in order to reduce the number of contour maps we have to compute the actual 
measure for. 

In order to calculate the measures, we need to integrate the joint density of the 
random field x. For example, P 2 is calculated as 

rb 1 rbn 

P r 2 — I ■ ■ ■ 7r(x) dx 

w cl\ J dfi 


where cq = Uj_ 1 and 5* = u e j +x , and for each i, j is the index that fulfills s* £ Gj. 
This high dimensional integral is generally difficult to calculate efficiently. How¬ 


ever, it can be estimated efficiently using the methods described in Bolin and 


Lindgren |2015j for latent Gaussian models with Markov properties, using sequen¬ 
tial importance sampling. For latent Gaussian models, the posterior distribution 
is generally not Gaussian, but the integral can still be estimated efficiently by us¬ 
ing, for example, the quantile correction method described in Bolin and Lindgren 

I 2 UI 5 I - 

In order to compute the P 0 measure, we also need to calculate the contour map 
function F u (-). Because the contour map function is an excursion function, we can 
again directly use the sequential integration method by Bolin and Lindgren 12015 


when calculating the contour map function. In order to use the sequential method, 
we need a parametric family for the possible contour avoiding sets for u. A natural 
choice of such a parametric family is D(a ) = {s : p( s) > 1 — a}, where p( s) are the 
marginal probabilities defined in 0- Thus, for a given value of a £ [0,1], D(a) 
defines a subset of f2 which serves as a candidate for the contour-avoiding set. See 
Bolin and Lindgren 2015 for further technical details of the method. 


3.2 Discrete and continuous domain interpretations 


For a problem defined on a discrete spatial graph, the crossings of a level can 
be described as the set of graph edges that connect nodes with values below and 
above the level. This does not quite match the usual concept of a contour curve, 
which requires the function to be defined at all spatial locations. It is therefore 
important to consider the link between the discrete computational model and the 
continuous spatial domain when interpreting the discrete calculations. 

For models representing cell averages on lattices, or other spatial sub regions, 
the natural interpretation of the discrete contour map calculations is to define the 
contour curves as the cell edges that separate nodes in different Gk sets. This is 


equivalent to applying Definition 2.2 to a piecewise constant function over continu¬ 
ous space. However, because we want to be able to handle more smooth functions, 
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we instead consider situations where the random field part of a model can be dis¬ 
cretized with weights for piecewise linear local basis functions. Common contour 
plotting methods are based on variations of such linear interpolation, e.g. contour 
in R |R Core Team 2013 and Matlab |MATLAB , 2014 . In practical generalized 


linear models, these continuous functions are combined with potentially discontin¬ 
uous covariates. The uncertainty about the resulting smooth level crossings and 
level-crossing jumps can be handled with more careful treatment of the sub-cell 
localization of the contour curves. 


3.2.1 Construction of continuous uncertainty interpretations 


Given contour map calculations on discrete point locations, the contour map func¬ 
tion F u (s) for the continuous spatial domain can be approximated, by assuming 
smoothness of the random field away from the contours. The main idea is to inter¬ 
polate the discretely computed F u (s) probabilities, assuming monotonicity of the 
random field between the discrete computation locations. For models constructed 
as piecewise linear basis approximations, this assumption is automatically fulfilled. 

The computational models in Section [4] and [5] use piecewise linear basis rep¬ 
resentations of stochastic PDE on triangles |Lindgren et ah, 20lT). Given the 
values of F u (s) on the vertices of a triangulation, the contribution to the 1 — a 
contour for F u (s) from each triangle T = (si,S 2 ,S 3 ) depends on if the triangle 
corners belong to the same or different G k level sets. If they belong to different 
sets, then the continuous version of F u (s) must go to zero somewhere inside the 
triangle, since it must have a contour separating the different Gk • In fact, this 
could be the case for any triangle, but triangles with common vertex Gk are un¬ 
likely to have interior contours. Because the precise location of interior contours 
are unknown and difficult to estimate, a practical approach to interpolating F u ( s) 
while taking Gk into account is to first eliminate all triangles with vertices that 
do not belong to the same Gk set. The 1 — a contour sets for F u (s) are then 
found by interpolation within each triangle of the pruned domain. The contour 
sets are made up of line segments for 1 — a crossings of F u (s) on the interior of 
the triangles, and partial sections of edges on the domain boundary for values 
of F u (s) above 1 — a. Optimal triangle interpolation based on the joint proper¬ 
ties of the entire model is intractable, but there are several parsimonious options. 
For interpolation to a location s = WiSi + u> 2 S2 + w 3S3 in the triangle T, where 
{(u>i, W 2 , w 3 ); wi,w 2 ,w 3 > 0, i Wk = 1} are barycentric coordinates, the op¬ 
tions in the excursions package include 


Step: F u (s) = min{F u (si), F u (s 2 ), F u (s 3 )}, 

Linear: F u (s) = Yll=i WkF u (s k ), and 
Log: F u (s) = exp{^'^ =1 u; fc log[F u (s fc )]}. 

For the Step and Log method, zero-width needle sets are avoided in the set con¬ 
struction by first eliminating any triangles where F u (s) = 0 for any of the vertices. 

A practical approach for visualizing the interpolated function is to compute a 
triangle subdivision, splitting each triangle into four new ones in each step, with 
the chosen interpolation method used to generate the values on the new vertices. 


12 























Figure 4: The left panel shows a triangulation with coloured indicators for the 
levels sets G 0 (red squares), G\ (blue circles), G 2 (cyan triangles) marked for each 
vertex. The middle panel shows results of the discrete calculations of F u ( s) at the 
vertices. Interpolation of the values within each triangle according to the algorithm 
in Section 3.2.1 followed by thresholding yields the contour avoiding sets shown in 
dark grey. The right panel shows the full interpolated version of F u (s), with the 
contour-avoiding set outline superimposed. 


Then, linear interpolation on the subdivided triangles, can be used when plotting, 
as illustrated in Figure |4} 

The qualitative behavior of the interpolation methods can be assessed by con¬ 
sidering test cases consisting of only two spatial locations, which allows exact 
calculation of F u (s) for a single level u along the line segment between the two 
locations. We define the joint distribution for x( 0 ) and x(l) to be 


x(0) 

x(l) 


~ N 


ho 

hi 


O’ o PO- O 0i 

poWi a\ 


and define the linear basis representation for s G [0,1] as x(s ) = (1 —s)x(0) + sx(l). 
The true F u (s) function is compared with its interpolation approximations for 
different combinations of distribution parameters (po, hi, oi), oh, h) and thresholds 
u. The three test cases shown in Figure [5] are 

(a) (ho, hi; of); o"i, p) = (0, 1, 1, 1, 0.9) with u — -0.5, 

(b) (/i 0 , hi, o~o, oq, p) = (0, 2, 4, 1, 0.9) with u = —0.1, and 

(c) (ho, hi> °o, en, p) = (0, 2, 4, 1, 0) with u = -0.1. 

The interpolation method Step is always conservative, Log is conservative for large 
target probabilities, whereas Linear is non-conservative for the target probability 
0.9 in test case (b). When both Linear and Log are non-conservative, Log stays 
closer to the true F u (s) function. In this simple example, x(0) and x(l) could be 
thought of as two of the neighboring locations on a mesh in two dimensions, and 
no major differences arise if one were to consider 2D interpolation over a triangle 
in the mesh. 


3.2.2 Numerical assessment of the continuous uncertainty interpreta¬ 
tions 


Coverage probabilities of credible contour regions constructed by different meth¬ 
ods were assessed by French 12014), for models where the true fields were on the 
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Figure 5: The true F u (s) (—), and interpolated approximations Step (—), Linear 
(-•), and Log for three test cases involving a model discretized at s = 0 and 
s = 1. A typical target probability 0.9 is shown for reference. The three test cases 
(a), (b), and (c) are obtained assuming a Gaussian distribution for [m(0) m(l)] T , 
where the parameters are varied (see the text in Section 3.2.1 for details) to show 
typical behaviors of the interpolation methods. One can note that the method Step 
is always conservative, Log is conservative for large target probabilities, whereas 
Linear is non-conservative for the target probability in case (b). 


same spatial resolution as the discrete models. The methods from Bolin and Lind- 


gren |2015 for a single contour level, which were generalized to multiple levels 
in Section 2.1 were assessed by performing discrete calculations based on the 
pointwise values at the lattice cell centers. Unfortunately, the comparison was 


designed to assess only the alternative contour concept from French |2014|, and 


lead to severely underestimated coverage probabilities for credible contour regions 
based on Definition |2.2| For an accurate assessment, one must either evaluate 


the pointwise calculations with respect to pointwise level exceedances, or take the 
properties of a continuous domain linear basis representation into account, as done 
in Section 13.2.11 

The numerical tests in Table [T| show that the interpolation methods are con¬ 
servative or on target when the smoothness assumption is fulfilled, but differences 
appear when the true random held is defined on a higher spatial resolution, where 
segments of the true contours can appear and disappear in between the discrete 
calculation points. More specifically, to assess the resulting credible contour map 
coverage probabilities for a single level u = 0, we simulated 100 zero mean random 
fields each from four different models, with the true model defined at a matching 
and higher resolution, and the differentiability v of the random held equal to 1 
and 2 in the model in Section |4j Each held was observed at 500 uniformly random 
locations, with the ratio between the variances for the held and the observation 
noise was hxed to 9, and this was repeated 10 times for each true held realization. 
The spatial range was hxed to 3, with the analysis domain as a 10 x 10 square, and 
the triangulation nodes for the discretized model were placed in a 20 x 20 lattice. 
The higher resolution helds were dehned on a 200 x 200 lattice. 

The resulting Monte Carlo estimates of the coverage probabilities for the in¬ 
terpolation methods Step, Linear, and Log, are shown in Table [lj for a target 
probability of 0.90. For comparison, the contour-avoiding sets for the nodes of the 
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Resolution 

V 

Step 

Linear 

Log 

Pointwise 

Matching 

1 

0.932(0.026) 

0.895(0.032) 

0.898(0.032) 

0.904(0.031) 


2 

0.932(0.026) 

0.887(0.033) 

0.891(0.033) 

0.919(0.029) 

High 

1 

0.850(0.037) 

0.583(0.052)- 

0.641(0.050)- 

0.901(0.031) 


2 

0.990(0.010)+ 

0.959(0.021)+ 

0.968(0.018)+ 

0.986(0.012)+ 


Table 1: Estimated coverage probabilities and associated Monte Carlo standard 


deviations for the simulation study in Section 3.2.2 The target coverage proba¬ 


bility was 0.90, and estimates significantly below or above the target are marked 
with — and +, respectively. For the Matching cases, the true field and numerical 
evaluations were on the same spatial scale, and for the High cases, the true field 
was defined on a higher resolution. The coverage probabilities for Step, Linear and 
Log are with respect to credible regions for continuous contour curves, whereas the 
coverage for the pointwise evaluation refers to the contour-avoiding sets of the field 
evaluated at the node locations of the discrete lower resolution model. 


discretized model were also evaluated, listed as Pointwise. When the true field is 
defined on a scale matching the discrete calculations, all methods have coverage 
probabilities matching the target, regardless of the smoothness of the continuous 
limit of the model. When the true field is defined on the higher resolution and 
with v = 1, the Pointwise evaluations are still on target, but the interpolation 
methods have lower coverage, due to not fully taking the within-triangle variabil¬ 
ity into account, and Step is the least non-conservative method, followed by Log, 
as expected from the test cases in Figure [5] In contrast, when the true field has 
v = 2, all four methods have overinflated coverage probabilities. This is likely due 
to the linear basis functions providing a relatively poorer approximation of the 
smooth covariance function, with the variance at each node of the low dimensional 


representation being larger than the variance in the continuous limit [see Lindgren 
et al., 2011, Bolin and Lindgren, 2013|. This gives higher uncertainty about the 


locations of crossings, which results in larger credible regions with higher actual 
coverage probability than the target. To avoid this, the computational resolution 
needs to be increased to more closely match the structural scale of the true field 
realizations. Further studies are required to give more precise guidance, but a 
basic ad hoc rule is that the longest edge of any triangle should be at most one 
tenth of the spatial correlation length for models with u = 1 but a quarter to a 
half may be sufficient for v = 2. The reason for this is that the true contour curves 
are rough for models with v = 1 but smooth for models with v = 2. 


4 An example using simulated data 


In this section, we illustrate the contour map methods using an example with 
simulated data. Let x(s), s e [0,10] x [0,10], be a mean-zero Gaussian Matern 
field. Its covariance function is given by 


cm i) 


2 1 ~ u (f> 2 

(47r)5r(z/ + | ) k 2v 


(K||h||)^( K ||h||), 


(4) 
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Figure 6: A simulated Matern field (left), the kriging prediction of the held based 
on 500 measurements (middle), and the corresponding standard errors (right). 
The uncertainty pattern seen in the standard errors is determined by the spatial 
pattern of the observation locations, and the kriging point estimate is smoother 
than the true held. 


Standard maps Pretty maps 


K 

1 

2 

3 

4 

3 

3 

5 

10 

Spacing 

1.577 

0.526 

0.394 

0.315 

2 

1 

0.5 

0.2 

Po 

0.613 

0.440 

0.267 

0.148 

0.616 

0.616 

0.407 

0.018 

Pi 

1.000 

0.999 

0.998 

0.966 

1.000 

1.000 

0.999 

0.291 

P 2 

1.000 

0.962 

0.523 

0.042 

1.000 

0.999 

0.878 

0.000 


Table 2: The quality measures Pq, P\ and P 2 for the contour maps shown in 
Figure [7J Here, K denotes the number of contours used and Spacing refers to the 
spacing between the contours. If the P 2 measure should be at least 0.9, one should 
among the Pretty contour maps use K — 3, and for the Standard maps one should 
use K = 2. 


where v is a shape parameter, k a scale parameter, (j) 2 a variance parameter, K u 
is a modified Bessel function of the second kind of order v > 0, and || • || denotes 
the Euclidean spatial distance. 

12011 


We use the SPDE representation by Lindgren et al. 


of the held and 

generate a realization of the held with parameters v — (j) — k — 1. The held is 
observed under additive N(0,0.001 2 ) noise at 1000 locations chosen at random in 
the square. Given these measurements the parameters and the marginal posterior 

20091. The sim- 


distributions are estimated using the INLA method Rue et al 


ulated held, the kriging prediction, and the corresponding standard errors can be 
seen in Figure [6] 

The Standard contour maps, based on the kriging predictor, with K — 1,... ,4 
contour levels are shown in Figure [7j as well as the hrst four Pretty contour maps. 
Table [2] shows the Pq, P\ and P 2 measures for the eight contour maps. If the P 2 
measure should be at least 0.9, one should among the Pretty contour maps use 
K = 3, which has a level spacing of 1. The third Pretty contour map, which has a 
level spacing of 0.5, has a P 2 measure just below 0.9. For the Standard maps, one 
should use K = 2, which gives a level spacing that is similar to the third Pretty 
map. 

The contour map functions F u ( s) for the standard contour maps can be seen 
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Figure 7: Standard contour maps (top) and pretty contour maps (bottom) for 
the simulated data example. As expected from the results in Table [2j the second 
Standard contour map is similar to the third Pretty contour map. 



Figure 8: Contour map functions F u (s) for the four Standard contour maps shown 
in Figure [7| with the contour lines of the corresponding contour maps superimposed 
in black. Dark green indicates areas of confidence. The spatial average of F u (s) 
gives the overall quality measure Pq, which is shown in Table [2] 


in Figure [8} Recall that the contour-avoiding set C u a can be retrieved from F u (s) 
as the set of locations where F u (s) > 1 — a. Also recall that the contour-avoiding 
set is the largest union of sets M* a C so that, with probability 1 — a, the field 
satisfies Uk < x(s ) < Uk+i for s e M* a . The F 0 measures corresponding to the 
contour map functions shown in Figure [8] are shown in Table [2] 

5 An application to temperature visualization 

In this section, we illustrate the contour map methods by looking at the problem 
of estimating the mean summer temperature for the US. The data we use is avail¬ 
able at www.image.ucar.edu/Data/US.monthly.met and was created from the 
data archives of the National Climatic Data Center. The data set contains daily 
measurements of min and max temperatures at approximately 8000 stations in the 
United States with a temporal coverage from 1895 to 1997. As a simple example, 
we focus on the mean summer temperature (June - August) for the year 1997. For 
each station, the mean temperature is defined as the mean of the recorded max 
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Figure 9: Mean summer temperature measurements for 1997. The spatial obser¬ 
vation pattern is fairly uniform with the exception of the mountainous region in 
the west. 


and min temperatures, and the values are averaged over the summer months to 
obtain the mean summer temperature. The resulting data is shown in Figure [9j 

As a first model, we assume that the temperature measurements are observa¬ 
tions of the true temperature field, x(-), under mean-zero Gaussian measurement 
noise, y l = x(sj) + e*, where e* are iid N(0,cr 2 ) variables. The latent field is mod¬ 
eled as x(s) = /3o + £(s) where /3 0 is the mean temperature and £(s) is a mean-zero 
Gaussian Matern field as in the example in the previous section. We estimate 
the model parameters using INLA and compute the kriging predictor of x(-) on a 
regular grid in the region. 

We calculate contour maps with K = 1,..., 10 levels based on the kriging 
predictor for the estimated model. The P-measures for the contour maps are 
shown in the left panel of Figure 10, and one can see that we are only allowed to 
use 2 contours if we want a contour map with P 2 > 0.9. 

However, modeling the temperature as a stationary Matern field is too simplis¬ 
tic. Therefore, we add altitude as a covariate for the mean of the field. Thus, the 
second model has x(s) = (3 0 + a(s)j3i + £(s) where a(s) is the altitude covariate. 
The altitude is available for each measurement location, and altitude data for the 
grid to which we do predictions was obtained from the ETOPOl 1 Arc-Minute 

20091. 


Global Relief Model Amante and Eakins 


The kriging predictor and standard errors using the second model can be seen 
in Figure [TTJ We again calculate contour maps and quality measures, now for 
K = 1,..., 20 levels. The resulting quality measures can be seen in the right panel 
of Figure 10 For this model, we are allowed to draw at most eight contours if we 
want a contour map with P 2 > 0.9. Figure [12] shows the Standard contour map 
with 10 levels and the corresponding contour map function. For this contour map, 
we have P 2 = 0.958 and P 0 = 0.206. 
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Figure 10: Contour map quality measures for the temperature application using a 
model without covariates (left panel) and using a model with altitude as a covariate 
for the mean (right panel). Squares indicate the P\ measure, triangles the P 2 
measure, and circles the P 0 measure. As expected from the definition, the Pi 
measure admits approximately twice as many contour levels as P 2 for a given target 
probability. The Pq measure has a different type of behavior, and the associated 
spatial pattern in Figure [12] is perhaps more useful. 



Figure 11: Posterior mean (left) and posterior standard deviations (right) for the 
temperature application using a model with altitude as a covariate for the mean. 
The spatial posterior mean pattern is dominated by the covariate, and the standard 
deviation pattern matches the spatial observation density seen in Figure [9] 


6 Discussion 


Although contour maps are widely used, relatively little research has focused on 
quantifying their statistical properties. We have here defined three contour map 
quality measures inspired by Polfeldt fl99 9 that can be used to assess how appro¬ 
priate a contour map is a for a given problem. The Pq measure is based on the 
spatially interpretable contour map function F u (-) and can be used as a measure of 
the proportion of the domain in which the held lies safely between the contour lev¬ 
els. For choosing an appropriate number of contour lines, or an appropriate level 
spacing, P 2 appears to be the most useful measure, as it has a precise definition 
that is also practically interpretable, ft gives the probability that the level cross¬ 
ing of the in-between levels associated with each Gk all fall inside their respective 
Gk region, making (G i,... ,Gk ) a joint credible region collection for those level 
crossings, as illustrated in Figure [fright). The Pi measure is similar in spirit to 
P 2 , but with no obvious advantages. 
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Figure 12: A contour map with eight contours with P 2 = 0.958 (left) and the 
corresponding contour map function F u ( s) with spatial average P 0 = 0.206 (right). 
As expected, the confidence is highest in the large flat regions in the central area 
due to the clear latitude effect. 


An intuitively interpretable approach to choose the number of equally spaced 
levels is to find K such that P 2 is above some threshold. For a joint credibility of 
90%, say, choose the largest K or the smallest spacing such that P 2 > 0.9. For a 
more permissive choice, require P 2 > 0.5 instead. In our experience, P 2 transitions 
quickly from 1 to 0, so the different thresholds are unlikely to give very different 
results. 

The computational methods are constructed on a discrete space, and the results 
have to be interpreted appropriately. In particular, we briefly presented methods 
to be used for continuous spatial domains, based on interpolation of the contour 
map function F u (s). We also showed that if the computational analysis is carried 
out on the same spatial scale as the true process, the credible contour regions 
resulting from the interpolation methods have the desired coverage probability. 
Further, when the true held has a fine scale structure, the analysis has to be 
carried out on a high resolution in order to capture small contour curve segments. 
We provided simple ad hoc rules for choosing the computational scale, but further 
studies are required to give more precise guidance. It should also be noted that 


the computational scale is of importance not only for the contour maps, see [Bohn 
and Lindgren |2011| for a more detailed discussion of this issue 


Efficient computational methods for finding the required quantities use the 
methods by Bohn and Lindgren |2015]. They work for the class of latent Gaussian 
models and are especially suitable for latent Gaussian Markov random held models. 
All methods discussed here are implemented in the R package excursions. 
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